/*********************************************************************************
 * Software License Agreement (BSD License)
 * Copyright (c) 2012, Carlos D. Correa (info@ngraph.org)
 * All rights reserved.
 * 
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 * 
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 * 3. Neither the name of netzen nor the
 *    names of its contributors may be used to endorse or promote products
 *    derived from this software without specific prior written permission.
 * 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * DISCLAIMED. IN NO EVENT SHALL <COPYRIGHT HOLDER> BE LIABLE FOR ANY
 * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 *  LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 *********************************************************************************/
/***************************************************************************
 *
 * Eigenvector Centrality Solver
 * based on IETL iterative solvers
 *
 * Copyright (C) 2009 by Carlos D. Correa <correac@cs.ucdavis.edu>
 *
 **************************************************************************/

#ifdef USINGLAPACK
#include "timing.h"

#include "node.h"
#include "centrality.h"


#include <boost/numeric/ublas/symmetric.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <ietl/interface/ublas.h>
#include <ietl/power.h>
#include <ietl/vectorspace.h>
#include <ietl/rayleigh.h>
//#include <ietl/lanczos.h>
#include <ietl/iteration.h>
#include <boost/random.hpp>
#include <boost/limits.hpp>
//#include <ietl/solver.h>

typedef boost::numeric::ublas::matrix<double> Matrix; 
typedef boost::numeric::ublas::vector<double> Vector;
typedef ietl::vectorspace<Vector> Vecspace;
typedef boost::lagged_fibonacci607 Gen;

#include "lafnames.h"       
#include LA_GEN_MAT_DOUBLE_H
#include LA_VECTOR_DOUBLE_H 
#include LA_SYMM_MAT_DOUBLE_H
#include "blaspp.h"
#include LA_SOLVE_DOUBLE_H
#include LA_GENERATE_MAT_DOUBLE_H
#include LA_EXCEPTION_H
#include LA_UTIL_H
#include LA_SOLVE_DOUBLE_H 
#include <lasvd.h>

inline void getDeriv(double *G, double *Qi, int i, int n) {
  for(int k=0;k<n*n;k++) Qi[k] = 0;

  for(int j=0;j<n;j++) {
    if(G[j*n+i]==0) continue;
    float sum = 0;
    for(int k=0;k<n;k++) {
      sum+=G[j*n+k];
    }
    if(sum>0) {
      for(int k=0;k<n;k++) {
		  Qi[j*n+k] = (G[j*n+k]/(sum+1.0) - G[j*n+k]/sum);
      }
      Qi[j*n+i] = ((G[j*n+i]+1.0)/(sum+1.0) - G[j*n+i]/sum);
    }
  }
}

EigenvectorCentrality::EigenvectorCentrality(Controller *controller):AnalysisTool("Centrality","EigenCentrality",controller) {
    setPropertyBool("Derivatives", false);
    setPropertyStringList("Edge Weight", "Default", "edge");
    setPropertyString("Result", this->getName());
	setPropertyStringList("Edge Filter", "Multiplicity", "edge");  /// Temporary. Figure out a general way to apply filters
}

void EigenvectorCentrality::compute() {
  assert(controller);
  // TO DO: set up these two via UI
  string resultProp = getPropertyString("Result");
  string resultDerivs = resultProp + ".Derivative";
  string resultDerivsT = resultProp + ".Derivative.Transpose";
  string resultDerivsM = resultProp + ".Derivative.Mutual";
	string resultDerivsB = resultProp + ".Derivative.Balance";
	string resultDerivsV = resultProp + ".Derivative.Magnitude";
	string resultDerivsN = resultProp + ".Derivative.Negative";
  Graph *g = controller->getGraph("graph");
  bool derivs = this->getPropertyBool("Derivatives");

  compute(g->nodes, derivs);
  controller->addDataProperty("node", resultProp, PROPERTY_FLOAT);
  g->computeNodeMinMaxFloat(resultProp);
  if(derivs) {
    controller->addDataProperty("edge",resultDerivs, PROPERTY_FLOAT);
	controller->addDataProperty("edge",resultDerivsT, PROPERTY_FLOAT);
	controller->addDataProperty("edge",resultDerivsM, PROPERTY_FLOAT);
	controller->addDataProperty("edge",resultDerivsB, PROPERTY_FLOAT);
	controller->addDataProperty("edge",resultDerivsV, PROPERTY_FLOAT);
	  controller->addDataProperty("edge",resultDerivsN, PROPERTY_FLOAT);
    g->computeEdgeMinMaxFloat(resultDerivs);
	g->computeEdgeMinMaxFloat(resultDerivsT);
	g->computeEdgeMinMaxFloat(resultDerivsM);
	g->computeEdgeMinMaxFloat(resultDerivsB);
	g->computeEdgeMinMaxFloat(resultDerivsV);
	  g->computeEdgeMinMaxFloat(resultDerivsN);
  }
}

int EigenvectorCentrality::compute(vector<Node*> &nodes, bool needDerivs) {
  int n = nodes.size();
	printf("Creating matrix of size %d x %d\n", n, n);
  double *G = new double[n*n];
  double *data = new double[n*n];
  double *characteristic = new double[n*n];
  double *qi = new double[n*n];

  string resultProp = getPropertyString("Result");
  string resultDerivs = resultProp + ".Derivative";
  string property = getPropertyString("Edge Weight");
  string resultDerivsT = resultProp + ".Derivative.Transpose";
  string resultDerivsM = resultProp + ".Derivative.Mutual";
  string resultDerivsB = resultProp + ".Derivative.Balance";
  string resultDerivsV = resultProp + ".Derivative.Magnitude";
  string resultDerivsN = resultProp + ".Derivative.Negative";
  print("Saving result as [%s]\n", resultProp.c_str());
	
	string edgeFilter = getPropertyString("Edge Filter");

  if(listener) listener->setProgress(0);

  for(int k=0;k<n*n;k++) {G[k] = 0; data[k] = 0;}

  for(int j=0;j<n;j++) {
    for(unsigned int k=0;k<nodes[j]->edges.size();k++) {
      assert(nodes[j]->edges[k]);
      assert(nodes[j]->edges[k]->getNode1());
      int i = nodes[j]->edges[k]->getNode1()->getIndex();
      float weight = 1.0;
      if(nodes[j]->edges[k]->hasProperty(property)) {
		  weight = nodes[j]->edges[k]->getPropertyFloat(property);
      }
		if(nodes[j]->edges[k]->hasProperty(edgeFilter)) {
			Property *p = nodes[j]->edges[k]->getProperty(edgeFilter);
			float value;
			if(p->getType()==PROPERTY_BOOL) value = (nodes[j]->edges[k]->getPropertyBool(edgeFilter))? 1.0:0.0;
			else value = nodes[j]->edges[k]->getPropertyFloat(edgeFilter);
			weight = (value > 0)? weight: 0;
			//if(value > 0) print("Property is %f\n", value);
			//else print("Property is %f\n", value);
		}
		G[j*n + i] += 0.5*weight; // this was	G[j*n + i] += weight; 
		G[i*n + j] += 0.5*weight;
    }
  }

  for(int j=0;j<n;j++) {
    float sum = 0;
    for(int i=0;i<n;i++) {
      sum += G[j*n+i];
    }
    if(sum>0) {
      for(int i=0;i<n;i++) {
		  data[j*n+i] = G[j*n+i]/sum;
      }
    }
  }

  for(int j=0;j<n;j++) {
    for(int i=0;i<n;i++) {
      double identity = (i==j)? 1.0:0.0;
      characteristic[j*n+i] = identity - data[j*n+i];
    }
  }


  timestamp t1 = now();

  // Initialize Matrix
  Matrix A(n,n);
  for(int i=0;i<n;i++) {
    for(int j=0;j<n;j++) {
      A(i,j) = data[j*n+i];
    }
  }
  
	/*
	printf("----\n");
	for(int i=0;i<n;i++) {
		for(int j=0;j<n;j++) {
			A(i,j) = data[j*n+i];
			printf("%f ", A(i,j));
		}
		printf("\n");
	}

	for(int i=0;i<n;i++) {
   float sum = 0;
    for(int j=0;j<n;j++) {
      sum += A(j,i);
    }
    printf("sum column %d is %f\n", i, sum);
  } 
	 */

  /*
  // Absolute Error Tolerance
  double abs_tol = 5. * std::numeric_limits<double>::epsilon();  
  // Relative Error Tolerance
  double rel_tol = std::numeric_limits<double>::epsilon();
  // Maximum Iterations
  int max_iter = 10*n;
  // Create the iteration object
  ietl::basic_iteration<double> iter(max_iter, rel_tol, abs_tol);
  // Create Vectorspace      
  Vecspace vec(n);
  
  // Create Generator for the starting vector
  Gen mygen;
  
  // Create the solver
  Solver<Matrix, Vector> mysolver;
   
  // Calculate the eigenvalue and the eigenvector
  std::pair<double, Vector> result = ietl::rayleigh(A, mygen, mysolver, iter, vec);
  
  // Print out the obtained results
  std::cout.precision(20);
  std::cout << "Eigenvalue: "<< result.first << std::endl;
  std::cout << "Eigenvector: "<< result.second << std::endl;
  
  // calculate the error as the norm of (Av-theta*v)
  Vector error = ietl::new_vector(vec);
  ietl::mult(mat,result.second,error);
  error -= result.first*result.second;
  std::cout << "error: " << ietl::two_norm(error) << std::endl;
  */

  /*
  typedef ietl::vectorspace<Vector> Vecspace;
  typedef boost::lagged_fibonacci607 Gen;  

  Vecspace vec(n);
  Gen mygen;
  ietl::lanczos<Matrix,Vecspace> lanczos(A,vec);

  // Creation of an iteration object:  
  int max_iter = 200;  
  double rel_tol = 50*std::numeric_limits<float>::epsilon();
  double abs_tol = std::pow(std::numeric_limits<float>::epsilon(),2./3);  
  std::cout << "Computation of 1 highest converged eigenvalues\n\n";
  std::cout << "-----------------------------------\n\n";
  int n_highest_eigenval = 1;
  std::vector<double> eigen;
  std::vector<double> err;
  std::vector<int> multiplicity;  
  ietl::lanczos_iteration_nhighest<double> 
    iter(max_iter, n_highest_eigenval, rel_tol, abs_tol);

  printf("[EIGEN] Computing eigenvector Lanczos method max_iters = %d rel_tol = %g\n", max_iter, rel_tol);
  try{
    lanczos.calculate_eigenvalues(iter,mygen);
    eigen = lanczos.eigenvalues();
    err = lanczos.errors();
    multiplicity = lanczos.multiplicities();
    std::cout<<"number of iterations: "<<iter.iterations()<<"\n";
  }
  catch (std::runtime_error& e) {
    std::cout << e.what() << "\n";
  } 
  for(int i=0;i<eigen.size();i++) {
    printf("eigenvalue[%d] = %f\n", i, eigen[i]);
  }

  std::vector<double>::iterator start = eigen.begin();  
  std::vector<double>::iterator end = eigen.begin()+1;
  std::vector<Vector> eigenvectors; // for storing the eigen vectors. 
  ietl::Info<double> info; // (m1, m2, ma, eigenvalue, residualm, status).
  try {
    lanczos.eigenvectors(start,end,std::back_inserter(eigenvectors),info,mygen); 
  }
  catch (std::runtime_error& e) {
    std::cout << e.what() << "\n";
  } 
  printf("[EIGEN] DONE\n");

  std::cout << "Printing eigenvectors:\n\n"; 
  for(std::vector<Vector>::iterator it = eigenvectors.begin();it!=eigenvectors.end();it++){
    std::copy((it)->begin(),(it)->end(),std::ostream_iterator<double>(std::cout,"\n"));
    std::cout << "\n\n";
  }
  */

  
  ietl::vectorspace<Vector> vec(n);
  boost::lagged_fibonacci607 gen;  
  int max_iter = 200; // orig 200
  double rel_tol = 1.0*std::numeric_limits<float>::epsilon();
  double abs_tol = std::numeric_limits<float>::epsilon();
//	double rel_tol = std::numeric_limits<double>::epsilon();
//	double abs_tol = std::numeric_limits<double>::epsilon();
	LaVectorDouble X(n);
  ietl::basic_iteration<double> iter(max_iter, rel_tol, abs_tol);
  
	
  printf("[EIGEN] Computing eigenvector Power method max_iters = %d rel_tol = %g\n", max_iter, rel_tol);
  std::pair<double,Vector> result = ietl::power(A, gen, iter, vec); 

  //std::cout << "Eigenvalue: "<< result.first << std::endl;
  //std::cout << "Eigenvector: " << result.second << std::endl;  

  Vector error = ietl::new_vector(vec);
  ietl::mult(A,result.second,error);
  error -= result.first*result.second;
  std::cout << "error: " << ietl::two_norm(error) << std::endl;

  printf("[EIGEN] Eigenvalue = %f\n", result.first);
  for(int i=0;i<n;i++) {
    X(i) = result.second[i]; 
    //printf("%g\n", X(i));
   }
  
	
	// Whoa! IETL failed with a simple 6x6 star network!!
	// Trying Lapack instead (but it may be up to 10x slower....
/*	
	LaGenMatDouble Mat(data, n,n);
	LaVectorDouble v(n);
	LaVectorDouble vimg(n);
	LaGenMatDouble Eigenvectors(n,n);
	
	for(int i=0;i<n;i++) {
		for(int j=i;j<n;j++) {
			Mat(i,j) = data[j*n+i];
		}
	} 
	
	LaEigSolve(Mat, v, vimg, Eigenvectors);
	{
		int i;
		int index = 0;
		double maxEig = -9999999;
		for (i = 0; i < n; i++) {
			if(v(i)>maxEig) {
				maxEig = v(i);
				index = i;
			}
		}
		
		i = index;
		
		printf ("eigenvalue[%d] = %g + %gi\n", i, v(i), vimg(i));
		for (int j = 0; j < n; ++j)
		{
			X(j) = fabs(Eigenvectors(j,i));
			//printf("%g\n", X(j));
		}
	}	
	 */
  /*
  LaGenMatDouble A(data, n,n);
  LaVectorDouble v(n);
  LaVectorDouble vimg(n);

  //LaGenerateMatDouble(A);
  LaGenMatDouble Eigenvectors(n,n);
  
  for(int i=0;i<n;i++) {
    for(int j=i;j<n;j++) {
      A(i,j) = data[i*n+j];
    }
  } 
  
  LaEigSolve(A, v, vimg, Eigenvectors);

  
 {
    int i, j;
    int index = 0;
    double maxEig = -9999999;
    for (i = 0; i < n; i++) {
      if(v(i)>maxEig) {
	maxEig = v(i);
	index = i;
      }
    }

    i = index;
    
    printf ("eigenvalue[%d] = %g + %gi\n", i, v(i), vimg(i));
    //printf ("eigenvector = \n");
    for (j = 0; j < n; ++j)
      {
	X(j) = fabs(Eigenvectors(j,i));
	//printf("%g\n", X(j));
      }
  }
  */
	printf("[EIGENVECTOR CENTRALITY] Time  %f sec\n", now()- t1);
  timestamp t2 = now();
  

  for(int i=0;i<n;i++) {
    nodes[i]->setPropertyFloat(resultProp,X(i));
    nodes[i]->createImplicitEdges(n, nodes);
  }

  if(needDerivs) {
    LaGenMatDouble Q(characteristic, n, n);
    LaGenMatDouble Qinv(n, n);
    LaGenMatDouble U(n,n);
    LaGenMatDouble VT(n,n);
    LaGenMatDouble D(n,n);
    LaGenMatDouble T(n,n);
    LaVectorDouble S(n);
    printf("[EIGEN] PseudoInverse of Q\n");
    LaSVD_IP(Q, S, U, VT);

    for(int i=0;i<n;i++) {
      for(int j=0;j<n;j++) {
		  D(i,j) = 0;
      }
      D(i,i) = fabs(S(i))>1e-6? 1/S(i): 0;
    }
    Blas_Mat_Trans_Mat_Mult(VT, D, T);
    Blas_Mat_Mat_Trans_Mult(T, U, Qinv);
   
     for(int i=0;i<n;i++) {
       if(listener) listener->setProgress((float) i/(float) n);
      printf("[EIGEN] Computing derivatives %d(%d)\r", i, n);
      getDeriv(G, qi, i, n);
      LaGenMatDouble Qi(qi, n, n);
      
      // text passed Qi is OK
		 /*
		 printf("Qi [%d]\n", i);
      for(int j=0;j<n;j++) {
        double sumQi = 0;
		for(int k=0;k<n;k++) {
		  sumQi+=Qi(k,j);
		  //printf("Qi[%d,%d] = %f\t%f ", j, k, Qi(k,j), G[j*n+k]);
			printf("%f ", Qi(k,j));
		}
		  printf("\n");
	  	//printf("  sumQi = %f\n", sumQi);
      }
		 */
      
      LaVectorDouble QiX(n);
      LaVectorDouble Di(n);
      Blas_Mat_Vec_Mult(Qi, X, QiX, -1.0);
      Blas_Mat_Vec_Mult(Qinv, QiX, Di, -1.0);
		 
		 
		 //printf("X = "); for(int k=0;k<n;k++) printf("%f ", X(k)); printf("\n");
		 //printf("QiX = "); for(int k=0;k<n;k++) printf("%f ", QiX(k)); printf("\n");
		 //printf("Di = "); for(int k=0;k<n;k++) printf("%f ", Di(k)); printf("\n");

      for(int j=0;j<n;j++) {
		  nodes[i]->getImplicitEdge(j)->setPropertySigned(resultDerivs, Di(j));
		  if(nodes[i]->hasEdge(j)) nodes[i]->getEdge(j)->setPropertySigned(resultDerivs, Di(j));

		  nodes[j]->getImplicitEdge(i)->setPropertySigned(resultDerivsT, Di(j));
		  if(nodes[j]->hasEdge(i)) nodes[j]->getEdge(i)->setPropertySigned(resultDerivsT, Di(j));

		  nodes[i]->getImplicitEdge(j)->setPropertySigned(resultDerivsV, fabs(Di(j)));
		  if(nodes[i]->hasEdge(j)) nodes[i]->getEdge(j)->setPropertySigned(resultDerivsV, fabs(Di(j)));
	
		  nodes[i]->getImplicitEdge(j)->setPropertySigned(resultDerivsN, -Di(j));
		  if(nodes[i]->hasEdge(j)) nodes[i]->getEdge(j)->setPropertySigned(resultDerivsN, -Di(j));
		  
	  }
	 }
	  printf("[EIGEN] Computing derivatives DONE\n");
      
	  timestamp t3 = now();
	  
	for(int i=0;i<n;i++) {
		for(int j=i;j<n;j++) {
			float di, dj;
			di = nodes[i]->getImplicitEdge(j)->getPropertyFloat(resultDerivs);
			dj = nodes[j]->getImplicitEdge(i)->getPropertyFloat(resultDerivs);
			float mutual = fabs(di)>fabs(dj)? di: dj; 
			float balance = fabs(di-dj);

			nodes[i]->getImplicitEdge(j)->setPropertySigned(resultDerivsM, mutual);
			nodes[j]->getImplicitEdge(i)->setPropertySigned(resultDerivsM, mutual);

			nodes[i]->getImplicitEdge(j)->setPropertySigned(resultDerivsB, balance);
			nodes[j]->getImplicitEdge(i)->setPropertySigned(resultDerivsB, balance);

			if(nodes[i]->hasEdge(j)) nodes[i]->getEdge(j)->setPropertySigned(resultDerivsM, mutual);
			if(nodes[j]->hasEdge(i)) nodes[j]->getEdge(i)->setPropertySigned(resultDerivsM, mutual);

			if(nodes[i]->hasEdge(j)) nodes[i]->getEdge(j)->setPropertySigned(resultDerivsB, balance);
			if(nodes[j]->hasEdge(i)) nodes[j]->getEdge(i)->setPropertySigned(resultDerivsB, balance);

		}
	}
	  
	  printf("[EIGENVECTOR DERIVS] Time  %f sec\n", t3 - t2);
  }
  
	
  
 
  delete G;
  delete data;
  delete characteristic;
  delete qi;
  if(listener) listener->done();

  return 1;
}



#endif //USINGLAPACK